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Optical carrier wave shocking: detection and dispersion 
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Carrier wave shocking is studied using the Pseudo-Spectral Spatial Domain (PSSD) technique. 
We describe the shock detection diagnostics necessary for this numerical study, and verify them 
against theoretical shocking predictions for the dispersionless case. These predictions show Carrier 
Envelope Phase (CEP) and pulse bandwidth sensitivity in the single-cycle regime. The flexible 
dispersion management offered by PSSD enables us to independently control the linear and nonlinear 
dispersion. Customized dispersion profiles allow us to analyze the development of both carrier self- 
steepening and shocks. The results exhibit a marked asymmetry between normal and anomalous 
dispersion, both in the limits of the shocking regime and in the (near) shocked pulse waveforms. 
Combining these insights, we offer some suggestions on how carrier shocking (or at least extreme 
self-steepening) might be realised experimentally. 



I. INTRODUCTION 

The self-steepening of an optical pulse envelope was 
first studied by DeMartini et al. in 1967 pL and is a 
well-known phenomenon associated with self-phase mod- 
ulation (SPM). Surprisingly however, the possibility of 
self-steepening of the optical carrier wave 1 was considered 
even earlier in a 1965 paper by Rosen ||], who showed 
that, for a third order x nonlinearity (and under suit- 
able conditions), a shock (or field discontinuity) develops 
in a finite distance. This latter phenomenon received lit- 
tle attention for more than 30 years, until it was revisited 
in the 1990s by Moloney et al. Q|, who performed Fi- 
nite Difference Time Domain (FDTD) simulations of the 
process. Three dimensional FDTD simulations of carrier 
shocking have recently been performed by Trillo et al. 

In the present paper, we investigate carrier wave shock- 
ing in nonlinear materials where SPM is accompa- 
nied by the generation of (odd) higher harmonics. For 
the dispersionless case we have generalized earlier pre- 
dictions based on the method of characteristics (MOC) 
to allow for arbitrary initial waveforms. This allows us, 
for example, to predict the carrier envelope phase (CEP) 
sensitivity of the shocking distance for optical pulses, as 
well as the dependence on pulse length. 

Our primary interest in this paper is the effect of linear 
dispersion on carrier shock formation. Since no analytic 
solutions exist in this case, we are forced to rely on nu- 
merical simulations. In any discussion of carrier wave 
shocking, it is important to distinguish between three 
related concepts: the physical system, the mathemati- 
cal model, and the numerical model. The mathematical 
model is an approximation to the physical system, while 
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1 Note that in this paper wc use the term 



carrier" to denote 
all the oscillations of the field, and do not use its other sense, 
i.e. that of fixed-frequency oscillations as used in a envelope and 
carrier representation of a pulse. 



the numerical model is an approximation to the mathe- 
matical model. Actual discontinuities occur only in the 
mathematical model, and it is these to which the idea of 
shock formation refers. In the physical system, disconti- 
nuities are prevented by phenomena not included in the 
mathematics, while numerical codes inevitably fail as a 
mathematical discontinuity is approached. The indicator 
of an imminent "shock" is the rapidly increasing gradi- 
ent (steepening) of the optical carrier. In the numerical 
model, we see numerical symptoms generated by extreme 
self-steepening, and these correlate with the onset of the 
discontinuities in the mathematical model. Under these 
circumstances, it has been necessary for us to develop 
a quantitative numerical test of "shock formation" . We 
have found the most satisfactory diagnostic to be "Local 
Discontinuity Detection" (LDD), which gives results that 
are in good agreement with theoretical MOC predictions 
based on the mathematical model. LDD provides a clear 
numerical measure of the rapid steepening that precedes 
the appearance of a discontinuity in the mathematical 
representation. We continue to use the term "shock" to 
describe the situation where a mathematical shock is im- 
minent. 

In our simulations, we exploit the flexibility in disper- 
sion management offered by the Pseudo-Spectral Spatial 
Domain (PSSD) technique to study carrier shock for- 
mation for a range of simple dispersion profiles, and de- 
termine the degree of phase mismatch that carrier shock- 
ing can tolerate. As expected, we find that shocking oc- 
curs when the nonlinearity dominates the linear disper- 
sion. However, it emerges that the process is asymmet- 
ric, with anomalous dispersion being far more conducive 
to shocking than normal dispersion. Hence, the fact that 
the LDD scheme does not detect carrier shocks in simula- 
tions involving (normally dispersive) fused silica, even at 
powers equal to its damage threshold, is neither surpris- 
ing nor necessarily discouraging. After all, anomalously 
dispersive materials could potentially be engineered. Fur- 
ther, our results also relate to how one might perform 
carrier shaping (as opposed to carrier steepening) , a pro- 
cess that has some interesting applications. 
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After briefly describing our simulation methods in 
section 0, we consider dispersionless shocking and the 
method of characteristics in and our LDD shock de- 
tection scheme in Then we discuss the effect of disper- 
sion on carrier shocking in followed by our numerical 
results in VI . In |VIl| we consider the potential relevance 



to experimental d etecti on of carrier steepening and/or 
shocks. Finally, in VII], we present our conclusions. 



II. SIMULATION METHODS 

The PSSD method [Q ||] offers significant advantages 
over the traditional FDTD and Pseudospectral Time- 
Domain (PSTD) H techniques for modeling the propa- 
gation and interaction of few-cycle pulses. Run times are 
generally faster, and PSSD also offers far greater flexi- 
bility in the handling of dispersion. Whereas FDTD and 
PSTD § propagate fields E{z),H(z) forward in time, 
PSSD propagates fields E(t),H(t) forward in space. It 
is important to keep this difference in mind when com- 
paring our results to those in ||. Under PSSD, the en- 
tire time-history (and therefore frequency content) of the 
pulse is known at any point in space, so arbitrary disper- 
sion incurs no extra computational penalty. In contrast, 
the FDTD or PSTD approaches use convolutions incor- 
porating time-response models for dispersion. 

We apply the PSSD algorithm to two representations 
of the field and source-free Maxwell's equations in non- 
magnetic media; the first uses the E and H fields, and 
the second the directional fields G ± (t) = a r (t) * E x (t) ± 
(3 r H y (t) Here the a r , (3 r include the (linear) per- 

mittivity and permeability of the material (i.e. e(i),/x). 
These G fields enable us to rewrite Maxwell's equations, 
and efficiently separate out the relevant forward-going 
part of the field. 

For an instantaneous x^ nonlinearity, the equations 
for E and H in the ID (plane wave) limit are 



dHy (t; z) 

dz 
dEJfi z) 

dz 



d_ 

' dt 



e Q e r (t)*E x (t-z) + x (Z) E x {t;zf 



= --^h*oHy(t;z)} 



a) 

(2) 



The G field simulations usually assume G~ = 0, and as 
a result contain only forward traveling components. The 
forward-only wave equation for G + is 



dG+(t;z) 
dz 



d 
dl 



f3 r a r (t) * G+(t; z) + l3 r x (3) E{t; z)\h 



where it is most straightforward to calculate the non- 
linear term by reconstructing E(t;z) from G + (t;z) in 
the frequency domain using E(u; z) = G + (uj; z)/2a r (uj), 
since G~ = 0. Notice the similarity between eqns. (pi) 
and (||), but that eqn. (||) propagates the field in a single 
first order equation, rather than two. 

Typical array sizes used in pulse simulations were 
N = 2 14 covering a time window T = 200fs, and (spatial) 



propagation steps were dz — OAcT/N « 0.9nm. We en- 
sured the stability of our integration using Orszag's 2/3 
rule pH , which involves setting the upper part of the 
spectral range to zero. It is worth noting that changing 
this cut-off, either by adjusting its position, or by using 
a smoothed (rather than step-like) filter, made little dif- 
ference to test simulations. The pulse profile used as an 
initial condition was 

E(t) = E sin(uj 1 t + 4>) sech(0.28wiVT), (4) 

where our standard parameters were lo\ = 2.356 x 
10 15 rad/s (i.e. A = 800nm) with r = 0.93. Such pulses 
are rather short (since the number of cycles inside the 
intensity FWHM is r) , but in fact the shocking distance 
is only weakly dependent on the pulse width, with signif- 
icant variation only appearing for pulses of a few cycles 
or less. 

We also performed CW simulations, for which we mod- 
eled just a single cycle of the carrier, assisted by the pe- 
riodic nature of the discrete Fourier transform. For these 
we used array sizes of N = 2 10 , and the time window was 
set by the period of the field oscillations. 

Our default value of nonlinear strength was x^Eq = 
0.02, which is comparable to that in fused silica at an in- 
tensity of 0.7 x 10 14 W/cm 2 ; our x^ parameter is equiv- 
alent to r\ in Rosen (2) , and a in Gilles ct al. Ji| . We use an 
instantaneous nonlinearity, since our primary interest is 
linear dispersion, and that is a far more significant effect 
than the nonlinear time response. 



III. CARRIER WAVE SHOCKING 

As an introduction to the process of carrier wave shock- 
ing, fig. ^ shows the profile of a pulse propagating in a 
dispersionless x medium, just before a shock occurs. 
The nonlinearity gives rise to a nonlinear index of refrac- 
tion U2E 2 , the effect of which is to increase the effective 
refractive index in the more intense regions of the pro- 
file. This reduces the phase velocity at the peak of each 
oscillation with respect to the rest of the waveform, and 
causes the slope on the trailing edges to increase dramat- 
ically. 

The effects seen in fig. [l] are associated with the gen- 
eration of third and higher harmonics, although the har- 
monic components are not particularly strong even when 
a shock is about to occur. As an example, fig. || shows 
how the harmonics build up as shocking is approached. 
The 10 scaling of the intensity spectrum in the figure ex- 
aggerates the contribution of the higher orders, and has 
been chosen for illustrative purposes. Notice that the 
profile becomes nearly flat (i.e. the spectrum falls off as 
the 4th power of the frequency) just before shocking is 
registered at around 4.3/xm. 

If we write 



E 



A(t) [sin(ojit) + 7 cos(3wit + ip)] , 



(5) 
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FIG. 1: The profile of a few-cycle optical pulse just prior to 
shocking in the dispersionless limit. The larger oscillations 
in the centre of the pulse undergo more self-steepening than 
those in the wings. The standard pulse parameters were used. 




FIG. 2: Development of the heights of the scaled harmonic 
peaks as a pulse approaches the (LDD) shocking distance of 
4.3^im in the dispersionless case. Each line perpendicular to 
the u/uii axis corresponds to the (scaled) contribution from 
that spectral peak. Note that the viewpoint has been rotated 
so that the contribution from the fundamental is to the right. 
The initial pulse contained about 33 cycles (r = 33). 



we find that an appropriate choice of A(t), with 7 = 
0.1 and tp = Q, gives us a passable match to fig. [|. 
Note that this choice of tp corresponds to the phase of 
third harmonic generation (THG) under index matched 
conditions, i.e. where no = n(u)\) = n(u}^). 

Rosen's original paper Q used the MOC to predict 
the formation of a value discontinuity in the field at cer- 
tain points within the profile. If the displacement of the 
dispersionless medium is written 



D = en 



E + x (1) E + x (3) E 3 



he showed that the wave equation for E is 



, d 2 E 
dx 2 



= 1+X 



(i) 



d 2 E 
dt 2 



■ X 



(3) 



d 2 E 3 

dt 2 ' 



(6) 



(7) 



and that the associated equation governing the charac- 
teristic lines of E is 



dE 
~dt 



v(E) 



dE 
dx 



(8) 



Here, the velocity v(E) is given by 
v(E) 



{e r + 3x^E 2 y/ 2 ' 



(9) 



where e r = 1 + x'' 1 '' = n a 2 is the (relative) dielectric 
constant and no the linear refractive index. 



v - dv 




FIG. 3: Method of Characteristics. Two points A and B 
on the field profile, separated initially by a time difference dt 
travel at different speeds v — dv and v, and meet at point C. 

Using eqn. (||) along with the construction shown in 
fig. [| we can derive a simple formula for the distance to 
shocking. The figure shows two characteristics AC and 
BC, originating from points A and B, and converging 
towards a shock at C after a distance of L. The intensity 
at A is higher than at B, so the speed associated with 
AC (represented by its gradient) is lower than that of 
BC. This means that at C, the field has two values, and 
a discontinuity has formed. From the geometry of the 
figure, it is easy to show that 



dv 



v 

T 



(c/n ) 



(10) 



where t, v = c/no, and L = vt are respectively time, 
speed and distance. 

On the other hand, differentiating eqn. (|9|) leads to 



dv 

a 



3cx 



(3) 



d(E 2 ) 



2 (n 2 + 3 X (3) E 2 ) 3/2 dt 

and this combined with eqn. ( JlCj ) yields 

cn Qy /l + 3x^E 2 /n 2 
3 X <- 3 H-dE 2 /dt) 
c y/l + 8n 2 E 2 /n 
4^ {-dE 2 /dt) ' 



(11) 



(12) 
(13) 



where n>2 — 3x /& n o is the material parameter deter- 
mining the intensity induced refractive index shift as 
niE 2 . For a given profile, a shock will occur first at 
the point where —dE 2 /dt reaches its negative extremum. 
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We can therefore define the shocking distance as 



L 



shock 



Min 



^l + 8n 2 E 2 /n 



c 
c 

Min-, — 

4n 2 (-dE 2 /dt) 



{-dE 2 /dt) 
1 



(14) 



for 8n 2 E 2 /n Q < 1(15) 



This formula is more general than that of either Rosen 
Pj or Gilles et al. H. Notice in particular that the 
parameter that controls the shock behaviour is not the 
gradient of the field (dE/dt), but that of the field squared 
{dE 2 /dt). 

For a sinusoidal initial waveform E = Eq sin(a;ii + 0), 
it is easy to show that shocks form at 



-it/4 + j-K — (f>, 



(16) 



where j is an integer, and that the shocking distance is 
2cuq c 



J shock 



3ujiX (3) E 2 ^m 2 E 2 ' 



(17) 



The approach is readily extended to pulsed waveforms 
of the type defined in eqn. (|J). Analytical results can 
be derived in the new situation in which eqn. (|l6|) be- 
comes a transcendental equation. However, the results 
are cumbersome and it is simpler to scan the profile nu- 
merically to determine the shocking parameters. It turns 
out that the shocking distance for short pulses exhibits 
an interesting sensitivity to the carrier envelope phase 4>. 
Whilst in the CW case all locations defined by eqn. JL6| ) 
were equivalent, for pulses, the one nearest the peak of 
the envelope has a shorter L s h oc k than the others. 

A set of results is displayed in fig. || where the shock- 
ing distance is plotted as a function of the carrier phase 
(j) for sech profiles with different pulse widths r. The 
dotted line is for the case of a very broad pulse for which 
L s hock is given by eqn. ( |l7| ). The sharp peaks mark a 
curve crossings where the shock location switches from 
one point on the sinusoid to another. Notice that the 
range of shocking distances increases for shorter pulses, 
and that, unsurprisingly, the lowest values occur when 
is around 7r/2 i.e. when the pulse has a cosine form. 

One very important point to note from fig. [| is that for 
pulses containing more than a few cycles, the dependence 
of the shocking distance on pulse width is very weak. We 
will see this message repeated later in section [v|, with 
shock regions being similar for both single cycle (r = 1) 
pulses and CW (r = oo) fields. 



IV. SHOCK DETECTION 

Since optical shock formation is directly associated 
with regions of increasingly steep field gradient, any nu- 
merical scheme, however sophisticated, is bound to fail at 
some point in the process. We therefore want to recognize 
when a shock is imminent, not only to avoid numerical 
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FIG. 4: MOC shocking distances as a function of phase, 
for pulses as in eqn. (^|), allowing for different pulse lengths. 
We include the 1/2 cycle r = 0.56 results to emphasize the 
trend; but even for a long pulse (r = 28), a small peak can 
be seen just below <f>/n — 0.01. The peak position is weakly 
r dependent. The + signs denote LDD shocking distances 
obtained from simulations of the r = 0.93 case. 



problems, but as a means to estimate the distance at 
which a discontinuity would occur in the mathematical 
model. 

One obvious symptom of impending numerical failure 
is loss of energy conservation |7|]. However, this does not 
give an indication at the first instance of a shock form- 
ing, but rather signals the accumulated effect of multiple 
small numerical failures from many shocked regions. 

A more physical strategy is to search for regions where 
the field gradient dE/dt is large and increasing rapidly, 
and to use this to predict the shocking point. A useful 
variant, suggested by the MOC calculation in the previ- 
ous section, is to use the value of —dE 2 /dt instead of the 
gradient. 

Overall, we find that the best method is Local Discon- 
tinuity Detection (LDD), which is similar to techniques 
used in other fields (see e.g. Jl^]). As the shock regime is 
approached, narrow shoulders with associated points of 
inflection appear within the regions of rapidly increasing 
gradient. The procedure is therefore to scan the field pro- 
file for the maximum gradient (of either E or E 2 ) and, if 
it occurs near a point of inflection, an incipient shock is 
registered. An example of a pulse that has just triggered 
the LDD diagnostic is shown in fig |[ 

The LDD method requires two parameters. The first 
determines the time scale used in determining whether a 
point of inflection exists. For this, we pick the scale set 
by our temporal grid and insist that the field gradients 
calculated at three adjacent grid points have opposing 
signs: either up-down-up, or down-up-down. The sec- 
ond determines the maximum range allowed between the 
maximum gradient and the point of inflection, and our 
default value for this was 10 grid points. In our simula- 
tions, we see that the position of the first detected shock 
depends only weakly on this range. We can easily mini- 
mize the small sensitivity to these parameters by holding 
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them fixed throughout any given set of simulations. As a 
result, we have found LDD to be a sensitive and reliable 
method of shock detection. 
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FIG. 5: The profile of a few-cycle optical pulse at the point 
of shocking. The larger oscillations in the centre of the pulse 
undergo more self-steepening than those in the wings, (a) The 
whole pulse, with the LDD carrier shock circled (o). (b) An 
expanded view of the shock region, with its point of inflection 
(at ~ 9.3fs) being very close to the steepest gradient (at ~ 
9.1fs). The standard pulse parameters were used. 

Eqn. (|l7]) predicts that the shocking distance should 
increase linearly with the refractive index no- We use 
this to test the LDD diagnostic in fig. q, where the ana- 
lytical formula is compared with the results of numerical 
simulations. This figure shows close agreement between 
prediction (using eqn. (|l5|)) and simulation, where the 
approximation causes the MOC prediction to be reduced 
by less than 0.1 um at no = 1. We see similar agreement 
between simulations using the LDD method and the car- 
rier phase sensitivity shown in fig. ||. The presence of 
small systematic differences (as on e.g. fig. [| ^) can be 
easily understood, since the LDD diagnostic is (strictly 
speaking) a test of the numerics, and is not a direct test 
for the presence of a physical shock or mathematical dis- 
continuity. 



THE EFFECT OF DISPERSION ON 
SHOCKING 



The primary purpose of this paper is to understand 
the principles of carrier shock formation in the presence 
of dispersion. Some simple ideas about the role of dis- 
persion can be understood from eqn.(^J) using the insight 
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FIG. 6: The LDD shocking distance as a function of refractive 
index, comparing the approximate MOC prediction from eqn. 
@ to PSSD simulations (+) for pulses with r = 0.93. Other 
PSSD simulation results from independent codes give very 
similar results to those shown on the graph. 



from eqns.(lO) - (jl5|), i.e. that the rate of change of E 2 
is the critical factor in shock development, rather than 
that of E. If only the fundamental and third harmonic 
components are considered (as in eqn.(||)), the effective 
refractive indices are 

An NL a = n 1 + ?l2 (l + 3 7 2 )/i, (18) 
An NL . 3 = n 3 + n 2 (3 + 7 2 ) Ii. (19) 

Evidently, the relative phase velocity of the two waves is 
affected by both linear and nonlinear dispersion, so the 
phase ip in eqn. (|5|) will vary accordingly as the pulse prop- 
agates. Fig. |?], which shows how dE 2 /dt varies with 
suggests that shocking is likely to be exacerbated when ip 
is small and positive, but moderated when ip is negative. 
Broadly speaking, the former case will be promoted by 
anomalous dispersion and the latter by normal disper- 
sion; however, the process is clearly complicated, since it 
involves time-dependent phase shifts between the waves 
and the interplay of linear and nonlinear dispersion. Of 
course, since carrier shocking relies on the establishment 
and maintenance of specific phase relationships between 
a set of harmonics, it must be expected that strong dis- 
persion of either sign will disrupt the shock formation 
process. On the other hand, the simple argument that 
has been offered suggests that shocking may be tolerant 
to a degree of anomalous dispersion, but not to a simi- 
lar amount of normal dispersion. In general, therefore, a 
graph of shocking signature versus refractive index mis- 
match might be expected to exhibit a shock region where 
nonlinearity dominates dispersion (displaced in the direc- 
tion of anomalous dispersion), surrounded by a shock- 
free region where dispersion dominates the nonlinearity. 
Moreover, if a dominant coherence length Lc can be de- 
fined, it is reasonable to expect shocking to occur when 
this exceeds the characteristic SPM length (Lspm)- As 
we shall see, all these features are borne out by the nu- 
merical results. 
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FIG. 7: The effect of different time-lags between fundamental 
and third harmonic on the steepness of the pulse. The field 
profile is as in eqn. (^), with A(t) = Aq and 7 = 1/16. (a) 
Scaled maximum values of Mj = Max^E 3 /dt) as a function 
of third harmonic phase offset tp. We can see that the maxima 
of dE/dt and dE 2 /dt occur at different if), (b) Position on the 
pulse 8j — ujitj of the maxima plotted in the top frame. Note 
that the kinks in (a) and (b) occur for offsets that give the 
longest shocking distances, not the shortest. 



In our numerical simulations, we make extensive use 
of model dispersion curves, which enable refractive in- 
dex differences between harmonics to be freely controlled, 
and lead to uncomplicated boundaries between the shock 
and no-shock regimes. The dispersion profiles shown 
in fig. |8] duly contain either refractive index steps An 
at the midpoints between successive odd harmonics, or 
smooth gradients 6 that provide a similar net change. 
In the simplest option, there is a single step at 2u)\, 
in which case the dominant coherence length is clearly 
Lc = 7r/ I&3 — 3fci| = TTc/3u>i\ri3 — rii\; we have also 
tried a single step at 4wi, which has a shorter Lc- In 
all cases, when n increases with frequency, the situation 
corresponds broadly to normal dispersion, while decreas- 
ing n corresponds to anomalous dispersion. The step 
size (or gradient) is chosen to be comparable to values in 
fused silica at u)\ — 1.5rad/fs, where the refractive index 
differences between the lower harmonics are Ani^ re 0.06 
and Att.3,5 ss 0.12W. While we consider the case of fused 
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silica itself in Section VII, the results based on the model 



dispersion characteristics are invaluable for understand- 



FIG. 8: Types of refractive index profile used. The solid lines 
show a single refractive index step midway between lui and 
W3, so that the fundamental is always phase mismatched from 
its higher harmonics, although they remain perfectly phase 
matched with each other. The dashed lines show multiple re- 
fractive index steps, each of the same size, and always midway 
between subsequent harmonics. The dotted lines show a lin- 
ear refractive index gradient which gives the same mismatch 
between subsequent harmonics as the multi-stepped case. We 
do not show lines for the case of a single step at 4wi to avoid 
cluttering the figure. 



ing the essential principles of carrier shock formation in 
dispersive media. 



VI. RESULTS AND DISCUSSION 

We will now analyze our numerical simulations of car- 
rier shocking in the presence of dispersion on the basis of 
the principles discussed above. Results are included for 
both a CW wave and for a single-cycle pulse. The CW 
case gives slightly wider shocking regions, but the differ- 
ences are minor. This is because the effect of the pulse 
envelope on the field amplitudes and gradients of the cen- 
tral carrier oscillation is small, except when considering 
sub-cycle pulses. In the results we present, the energy 
conservation and LDD measures reveal the presence of 
sharp boundaries between the shocking and non-shocking 
regions. 

Results for single refractive index steps at 2wi and 4wi 
are presented in figs. ^| and |l^ respectively, whereas fig. 
[Tl| has a step midway between all harmonics. In both 
single stepped cases, a useful coherence length can be 
defined on the basis of the mismatch between the fun- 
damental and the harmonic just above the step. In the 
multi-stepped case, there is no easy way to define a dom- 
inant coherence length. 

Fig. ^(a) shows the pulse profiles in the anomalous 
(negative) and normal (positive) single step cases, for the 
smallest step size at which shocking did not occur. The 
first obvious difference between the profiles is the oppo- 
site walk-off direction of the third harmonic, although the 
third harmonic contribution is hard to see in the anoma- 
lous case. The second is that the pulse profiles exhibit 
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FIG. 9: Carrier shocking for CW and pulsed (r = 0.93) cases 
with a single refractive index step at 2wi. Frame (a) compares 
chirped pulses at a 40/im propagation distance from the non- 
shock region immediately outside of the shocked region; the 
upper curve (X) is for the negative step, the lower (Y) for 
the positive step. Frame (b) shows the correlation in the 
region where Lc < Lspm between energy conservation failure 
[logarithmic right hand scale, dots) and the LDD detected 
shocking distance (left hand scale, solid line). The dotted 
line shows the LDD results for a CW field (r = oo). 



distinctly different characteristics according to the sign of 
the dispersion. Narrow spikes are visible in the anoma- 
lous case, whereas profiles with more rounded maxima 
occur for normal dispersion. A possible interpretation is 
suggested by fig. [?], where we saw that anomalous dis- 
persion tends to create regions of higher gradient. 
The shocking region in fig. 



10 has a similar outline 
to that of fig. ||(b) except that it is slightly narrower, 
as expected from the coherence length discussion above. 
However, the shocking region in fig. |Il] is much nar- 
rower, especially given the change in the Art scale. This 
is because the multi-stepped nature of the refractive in- 
dex leads to a correspondingly large range of coherence 
lengths, with shorter ones corresponding to those span- 
ning several steps. Since the shocking region in the multi- 
step case has reduced in size by a factor of two or three, 
a reasonable inference might be that the dominant co- 
herence length results from interaction over two or three 
refractive index steps. 

In realistic media, the refractive index will vary 
smoothly with frequency, and the group velocity can dif- 
ferent from the phase velocity. We can approximate this 
situation most simply using a refractive index gradient 




FIG. 10: Carrier shocking for CW and pulsed (r = 0.93) 
cases with a single refractive index step at Alo\\ showing the 
correlation in the region where Lc < Lspm between energy 
conservation failure [logarithmic right hand scale, dots) and 
the LDD detected shocking distance (left hand scale, solid 
line). The dotted line shows the LDD results for a CW field 
(r = oo). 




-0TJ5" 
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FIG. 11: Carrier shocking for CW and pulsed (r = 0.93) cases 
with multiple refractive index steps; showing the correlation 
in the region where Lc < Lspm between energy conserva- 
tion failure (logarithmic right hand scale, dots) and the LDD 
detected shocking distance (left hand scale, solid line). Note 
the narrower range of An as compared to the previous two 
graphs. The dotted line shows the LDD results for a CW field 
(r = oo). 



rather than a series of steps. The results using the LDD 
method in this case can be seen on fig. |lj, where now 
we also vary the strength of the nonlinearity. As in the 
previous cases, we see a well defined shocking regime that 
is asymmetric about the non-dispersive case. 

Detailed examination shows that the curves in fig. [l^ 
exhibit a marked similarity, and can be brought into 
near perfect coincidence by applying the scaling x^ ~ > 
X^/m, L — > mL, and An — ► An/m. We also get a 
comparable similarity for each of figs. 0,nQ,O, when sim- 
ulations at nonlinear strengths of x Eq = 0.01 and 
0.04 are added to the results. This demonstrates that 
the character of the shocking is dominated by the sign of 
inter-harmonic phase velocity differences, not by the local 
group velocity dispersion at each harmonic. Thus anoma- 
lous (normal) dispersion is primarily interesting because 
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it gives a negative (positive) refractive index shift be- 
tween successive harmonics. We leave the complicated 
(and far more subtle) effects of group velocity differences 
or dispersion for later work. 
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FIG. 12: Shocking distance vs mismatch for weak refractive 
index gradients 8, and a range of nonlinearities. This shows an 
abrupt cut-off at positive An, but a relatively gradual one at 
negative An. The initial pulses are identical to those in fig. [)| 
No CW results are shown, since the dispersion experienced by 
the field would be identical to the multi-stepped case shown 
in fig. 0. 

An important feature of all the results is the pro- 
nounced asymmetry, with shocking persisting much fur- 
ther into the anomalous dispersion regime; we have even 
seen it for the case of a weakly parabolic refractive index. 

We can quantify this asymmetry by considering lin- 
ear and nonlinear contributions to the phase matching of 
the harmonics (i.e. linear and nonlinear refractive index 
shifts), as described at the beginning of section^. In the 
simple single step case shown on fig. ^|(b), the lst-3rd 
harmonic phase shift will dominate, because it applies to 
the two most intense spectral components. We can calcu- 
late the SPM-induced refractive index shift between the 
fundamental and third harmonic with eqns. @ and (U) 
to be An NL: i-An NL:3 = -0.011, since x (3) £^= 0.02 at 
the peak of the carrier oscillations, the refractive index is 
nl = 2, and 2n%E\ = 2 x 3x (3) E 2 /8n . Tni s is roughly 
comparable to the offset of the shocking region, which 
is centred at about Ani^ ~ —0.014. We cannot ex- 
pect perfect agreement, since the calculations ignore the 
role of higher harmonic generation and depletion of the 
fundamental. Whilst a simple calculation is reasonably 
successful in this single-step case, it cannot be applied 
for a realistic medium - or indeed to the situation shown 
in fig. [IJ. There, the effects of dispersion, higher har- 
monic generation, and nonlinear refractive index shifts 
are inextricably intertwined. 

To summarize, we have demonstrated that shocking 
is strongly dependent on the interplay between Lc and 
Lspm, with shorter Lc's (increasing An's) decreasing 
the likelihood of shocking. We have also deduced the 
reasons for the strong asymmetry of the shocking region. 



VII. APPLICATIONS AND EXPERIMENTS 

In sections [j] and [Tv| , we discussed how imminent car- 
rier shocking might be recognised computationally. In 
considering whether shocking might be detectable exper- 
imentally, we must now decide how it might manifest 
itself in the laboratory. A mathematical discontinuity is 
clearly not a physical possibility, even if Rosen @ did 
manage to accommodate it theoretically, albiet at the 
expense of energy conservation. In practice, the increas- 
ing field gradients (and spectral broadening) that precede 
a shock will inevitably engender new physical processes 
that will limit the steepening. Indeed, we have already 
seen this happening in the previous two sections, where 
dispersion has been seen to frustrate the self-steepening 
process; the next barrier would be the time-scale of the 
nonlinear response. 




FIG. 13: Numerically predicted shocking distances L s hock 
in fused silica. These are the MOC predicted distances for 
the waveform, assuming the dispersion was (abruptly) ne- 
glected, the shortest distance shown is about 1.9/im. The 
damage threshold for fused silica was taken to be Pthreshoid = 
50TW/cm 2 . 



Although there is no question of a mathematical dis- 
continuity being observed in an experiment, the recent 
advances in the measurment of optical pulse profiles (see 
e.g. p3[ ), suggest that it might be possible to observe 
carrier steepening. Indeed, it has already been pre- 
dicted that noticeable steepening effects could occur 
for pulses propagating in fused silica. Unfortunately, in 
our own simulations of this process, the LDD shock de- 
tection was not triggered at any realistic pulse intensity. 
Although visible steeping did occur, at no point was a 
seriously distorted waveform approached, and the effects 
would have been milder if we had included the finite re- 
sponse time of the nonlinearity. The results presented 
in fig. [l3] show the steepest gradient recorded as a func- 
tion of distance and pulse intensity for these simulations. 
The third harmonic coherence length for these parame- 
ters is about 7fim, and we can attribute the regular vari- 
ation with distance seen in the figure to the third har- 
monic component aligning with successive oscillations of 
the fundamental as the waves move across each other. 
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The key reason why incipient shocking was not detected 
is because the pulse frequency lies in a region of weak 
normal dispersion, whereas we have shown in section |v| 
that a region of weak anomalous dispersion (assuming it 
could be achieved) would be favourable. 

The current interest in media with tailored dispersion 
characteristics [jl4], |l5, [l6| raises the interesting question 
of whether it might be possible to engineer a material 
that maximize self-steepening. A major stumbling block 
would be the need for control over many harmonic orders; 
however, for a narrow-band pulse, only the dispersion 
characteristics close to the harmonics would be relevant, 
which might perhaps make the technical challenge less 
formidable. 

The recognition that dispersion control can enhance 
(or reduce) carrier self-steepening suggests other appli- 
cations if we widen our horizons to encompass the more 
general idea of carrier shaping. In this case, we would 
exploit both nonlinearity and dispersion control to opti- 
mize the shape of the carrier oscillations for a particular 
experiment. Applications such as high harmonic genera- 
tion (e.g. |l7j]) might well benefit from suitably designed 
carrier wave modulation. 



VIII. CONCLUSIONS 



merical simulations, and shown how it is influenced by 
a range of parameters, particularly dispersion. We have 
also obtained remarkable agreement between numerical 
simulations and theoretical predictions of the shocking 
distance in the dispersionless limit, and shown that the 
process is sensitive to both CEP and pulse duration. 

Although we have confirmed that shocking occurs in a 
narrow parameter range, this is far from being the whole 
picture. In particular, there is a distinct asymmetry be- 
tween the anomalous and normal dispersion regimes. The 
former leads to shocking signatures such as the appear- 
ance of narrow spikes as the higher harmonics interfere 
on the steepening part of the pulse profile. In contrast, 
normal dispersion creates no such features, and the pulse 
profiles have a rather blunt appearance. The asymmetry 
arises from the effect of the nonlinear refractive index on 
the dispersion induced phase mismatch. 

The conclusion to be drawn from our results is clear: if 
they could be engineered, materials with wide regions of 
weakly anomalous dispersion are much better candidates 
for generating steep, shock-like field profiles than those 
(such as silica) with a weak normal dispersion. Detecting 
incipient shock formation in materials like fused silica is 
likely to be a near impossible task, given the constraints 
imposed by their damage thresholds. 



In this paper we have investigated carrier shock for- 
mation, developed criteria for detecting its onset in nu- 
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